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Starting from a mean-field model of the Bose-Einstein condensate dimer, we reintroduce classi¬ 
cally forbidden tunneling through a Bohr-Sommerfeld quantization approach. We find closed-form 
approximations to the tunneling frequency more accurate than those previously obtained using dif¬ 
ferent techniques. We discuss the central role that tunneling in the self-trapped regime plays in 
a quantitatively accurate model of a dissipative dimer leaking atoms to the environment. Finally, 
we describe the prospects of experimental observation of tunneling in the self-trapped regime, both 
with and without dissipation. 


A Bose-Einstein condensate of atoms in two modes 
(BEC dimer) is a simple interacting quantum system 
that has recently become accessible to increasingly pre¬ 
cise experiments [Tj. It has been used to demonstrate 
matter-wave interferometry [5], number squeezing EMI 
and measurements transcending the standard quantum 
limit HE], and its prospective applications include grav¬ 
ity detectors [5j, noise thermometers [9J and tests of the 
EPR paradox m- Especially exciting is the opportunity 
to study the gradual emergence of classical mechanics as 
the number of particles in the system is increased HU. 

The simplest theoretical approach to the BEC dimer 
is the two mode mean-field model m nsj. In this 
model, a phenomenon known as self-trapping takes place: 
a coherent state prepared in certain regions of phase 
space remains in the neighborhood of the nearest stable 
fixed point (self-trapping point) forever. Self-trapping 
has been experimentally observed for relatively short 
times mm- In the quantum treatment of the two¬ 
mode model, tunneling between the two self-trapping 
fixed points eventually occurs. This process of “quan¬ 
tum sloshing” generates macroscopic entanglement be¬ 
tween the two wells of the dimer m- The time scale on 
which tunneling takes place can be found numerically by 
directly integrating the Schrodinger equation, but this of¬ 
fers little insight into the process. An analytical estimate 
of the tunneling frequency has been obtained using quan¬ 
tum perturbation theory (see m and references therein), 
but the expansion employed is only valid in a parame¬ 
ter range where the tunneling frequency is exponentially 
small. 

In this paper, we use the semiclassical quantization 
approach to the Bose-Hubbard dimer pioneered in m 
to obtain highly accurate analytical approximations to 
the tunneling frequency of the two lowest-energy states. 
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Unlike quantum perturbation theory, the semiclassical 
techniques remain applicable as long as approximately 
self-localized quantum states exist. This showcases the 
power of the semiclassical approach to many-body prob¬ 
lems, and allows us to clarify the dependence of the tun¬ 
neling time on the system’s parameters. We also discuss 
the prospects of experimentally observing tunneling in 
the self-trapped regime. 

The Bose-Hubbard dimer is mathematically equiva¬ 
lent to a spin system and to a certain limit of the 
Lipkin-Meshkov-Glick model. Semiclassical quantiza¬ 
tion of these equivalent models was considered by [18l - 
[M! and [21] . respectively. We complement these ear¬ 
lier works by providing a connection to the quantization 
condition of m, proposing closed-form approximations 
valid in the relevant parameter range and offering a dis¬ 
cussion of cold-atom experiments that could probe the 
tunneling phenomenon. 

The rest of the paper is organized as follows. In Sec¬ 
tion [I] we review the Bose-Hubbard dimer and its mean- 
field approximation. Section [TT] is devoted to tunneling 
between the fixed points using exact diagonalization re¬ 
sults. Iu Section III we introduce the semiclassical quan¬ 
tization condition and obtain a closed form expression for 
the tunneling frequency. Finally, in Section [TV] we discuss 
applications of this expression to problems of entangle¬ 
ment and atom loss rate from a dissipative optical lattice, 
as well as prospect of experimental confirmation. 


I. THE BOSE-HUBBARD DIMER 

We will consider bosonic atoms in a double well op¬ 
tical trap sufficiently deep that only the lowest state in 
each well is populated. In this so-called two mode ap¬ 
proximation, the atoms’ dynamics is described by the 
Bose-Hubbard Hamiltonian [22] , 

H = - J(a\a 2 +a\ai)+^ (ni(hi - 1) + h 2 (h 2 - 1)) (1) 
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where at is the annihilation operator for a boson in well 

1 and n t = a\a,i is the number operator. 

Of special interest are the coherent states of the 
model [23} (24j. These states correspond to all atoms be¬ 
ing a single BEC [25j [263 and can be characterized by 
their expectation values of the population imbalance be¬ 
tween the wells, z = (A/ — TVs)/2, and of their relative 
phase, (j). In terms of the creation operators, the coherent 
states can be expressed as, 

k, </>) = (j{l + z)/2a\ + V(l-z)/2e^4) N |0). 

( 2 ) 

For large numbers of particles, the coarse dynamics of 
these states is well approximated by a bosonic Josephson 
junction (BJJ) model in which z and <j> are the dynamical 
variables. The Hamiltonian of this model is, 

H = —^-\/l — z 2 cos (/>, (3) 

where A = and the dimensionless time is r = 

2 Jt/h [TB]. The BJJ model exhibits a bifurcation [27] 
at A = 1: as A is increased beyond this critical value, a 
stable center at z = 0, </> = ir breaks down into a sad¬ 
dle point point at the same coordinates and a pair of 

stable centers at z = ±^j 1 — , (j) = n. These stable 

centers, corresponding to a persistent population imbal¬ 
ance between the dimer’s two wells, are known as the 
self-trapping points. 


II. TUNNELING BETWEEN THE 
SELF-TRAPPING POINTS 

Within the BJJ model, the self-trapping fixed points 
are stable: a trajectory initially sufficiently close to one 
of them remains close to it for all time. In the full Bose- 
Hubbard dynamics, however, tunneling between the two 
self-trapping points occurs with a finite frequency. An ex¬ 
ample of this process is shown in Figure [I] which depicts 
the Husinri function 128] , a quasiprobability distribution 
over the coherent states \z, </>) given by, 

Qv-k, <t>) = Ik, <^#)| 2 ( 4 ) 

for a pure state \ip). The Husimi function is initially cen¬ 
tered at one of the fixed points, but over time it tunnels 
to the other, and then back again. 

A quantitative signature of the tunneling is an oscil¬ 
lation of the wells’ populations. The frequency of this 
oscillation can be found by numerically integrating the 
Schrodinger equation of the Bose-Hubbard dimer for a 
long time and computing the power spectrum of the well 
populations. The most prominent feature in the spec¬ 
trum corresponds to the tunneling frequency. 

Since the dynamics of the coherent state near the self¬ 
trapping fixed points appears very simple, we may try 


t = 0 0.52s 1.05s 1.57s 2.10s 



0.8 1 1.2 

t) I n 


FIG. 1. Tunneling between the self-trapping fixed points. 
In the BJJ model, trajectories sufficiently close to the self¬ 
trapping point remain confined to its neighborhood for¬ 
ever (far left panel). However, as shown in the remain¬ 
ing panels, in the Bose-Hubbard model the Husimi func¬ 
tion of a coherent state initially centered at the z > 0 self¬ 
trapping fixed point tunnels from one fixed point to the 
other. For the full video from which these stills are taken, 
see http://youtu.be/hX4nhoMb4G0, Parameters: N = 40 
atoms, with A = 1.1 and J = 10 Hz. 
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FIG. 2. The probability of observing the coherent state cen¬ 
tered at the self-trapping fixed point in one of the n most 
probable states, for n = 2, 3, 4,..., as a function of the par¬ 
ticle number N. (A = 1.025, U = 2n x 0.063 Hz.) 


to reduce the dimensionality of the problem by restrict¬ 
ing the system to some subspace of the Hilbert space. 
Remarkably, in the neighborhood of the mean-field fixed 
points, only a few energy eigenstates contribute appre¬ 
ciably to the coherent state eh uni- How many states 
need to be accounted for depends on the particle num¬ 
ber (see Figure [2j. Our intuition is that as N increases, 
the “size” of the coherent state in phase space shrinks, 
but the “size” of the eigenstates shrinks even faster, and 
ever-more eigenstates are needed to correctly account for 
the coherent state dynamics. However, even for a few 
hundred atoms much of the tunneling dynamics can be 
captured by keeping just two states (see Figure|3|. At the 
self-trapping fixed points, these two states are the pair 
of highest energy states of the Bose-Hubbard model Eli. 
They are symmetric and antisymmetric combinations of 
states localized in each well. 

The energy splitting between the symmetric and 
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FIG. 3. A two-state description of the tunneling remains 
valid as N increases, although the dynamics is more com¬ 
plex as the system becomes less discrete. The Husimi 
function is shown at the five times spaced by a quar¬ 
ter of the tunneling period expected from the two-state 
model. For the full video from which these stills are taken, 
see http://youtu.be/p_LL85VBohU, (N = 500 atoms, A = 
1.025 and U = 2n x 0.063 Hz.) 




FIG. 4. The BJJ result of zero tunneling frequency for A > 1 
is gradually approached by the two-eigenstate model as the 
number of atoms increases. Nonetheless, a nonzero frequency 
is expected for any A and any finite N. In all plots, J = 10 Hz. 


antisymmetric states agrees closely with the oscilla¬ 
tion frequency extracted by numerically integrating the 
Schrodinger equation [52] , The splitting between these 
states can also be computed for A < 1; in this case, there 
is only one fixed point at (f> = 7r, and the energy split¬ 
ting closely agrees with the BJJ frequency of oscillations 
about that point. Both above and below A = 1, the BJJ 
limit is approached as N is increased (see Figure [4]). 

The energies of the two highest-energy states are easily 
found numerically even for very large N, but it is desir¬ 
able to explain the simple trends with N and A shown 
in Figure [4] using an analytical model. Quantum per¬ 
turbation theory can be used to obtain estimates of the 
tunneling frequency for small J/U ~ N/A |T6j[301 1331134] , 
but not in the region A « 1 where tunneling becomes a 
significant effect. In the next section, we will pursue an 
alternative approach. 


III. SEMICLASSICAL QUANTIZATION 

To shed light on the convergence of the results of the 
two-state model to those of the BJJ, we will start with the 
BJJ model and recover additional features of the dynam¬ 
ics through Bohr-Sommerfeld quantization. Graefe and 
Korsch m applied Bohr-Sommerfeld quantization to 
this problem numerically, obtaining excellent estimates 
of the eigenenergies even for atom numbers N < 10. 
In this section, we start from their formulation of the 
quantization condition but proceed analytically to pro¬ 
duce accurate closed-form expressions for the tunneling 
frequency. 

The quantization condition in the self-trapping region 
of the symmetric dimer described by the Hamiltonian of 
Eq. []TJ: is pTj, 

Vl + k 2 cos(2 S w — St) = —k. (5) 

Here, 2 S w is the action associated with the self-trapped 
classical orbit, k = exp(—7r S e ) and 2 S e is the (Euclidean) 
action associated with tunneling. Both S w and S e are 
measured in units of Planck’s constant, h , and so are 
dimensionless [35]. The phase correction term St can be 
expressed in terms of S € as m 

St = arg r Q + iS^ - S e In \S e \ + S e . (6) 

For a discussion of the physical significance of S</,, see [SB] 
pp. 50-51]. 

The actions S w and S e are functions of the energy E 
and the nonlinearity A, and can be expressed as inte¬ 
grals over phase space (see Figure [5]); this is discussed in 
greater detail in Appendix [B] 

Let us assume that the energy splitting between sym¬ 
metric and antisymmetric combinations of states local¬ 
ized in the two self-trapping regions of phase space is 
small relative to the spacing of allowed energies in each 
region. As shown in Appendix [A] in this case the quan¬ 
tization condition implies the splitting is approximately 

A E= — exp(nS e ), (7) 

7 r 

where lu is the frequency of the classical motion in a 
self-trapped orbit (related to the action of the orbit S w , 
since 2ir/hu! = T/h= 2 dS w /dE) and S f is as before the 
Euclidean action associated with the tunneling. These 
quantities depend on the shape and size of the classical 
orbits, which are determined by A and the energy of the 
unperturbed state E. 

Let the classical turning points be z± (see Figure [b]). 
The size of the orbits is captured by the dimensionless 
parameter, 


k = 



( 8 ) 
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FIG. 6. Pairs of classical orbits and their turning points z±. 
The orbits on the left (A = 2, E = 1.15) are librations, while 
those on the right (A = 4, E = 1.15) are rotations. 


FIG. 5. The actions appearing in the quantization condition 
[Eq. 0] have a geometric interpretation. This figure depicts 
the phase space of the BJJ model for A = 2. The grey curves 
are trajectories; the actions S w and S t for energy E = —1.15 
are equal to the areas of the marked regions. In the case of 
S w , the action corresponds to the phase space area of the 
classical orbit. 


Furthermore, let, 

k' = \J\ — k 2 = —, and a 2 = X- ~. 

z+ z\- 1 

In Appendix [B] we show that in terms of these quantities 
the splitting A E of the highest-energy state is given by, 


A E = — exp(7r5' £ ) 

= exp (r (iV +1} (■ ( x - x) k ' ] + *+ (K(fc,) - E(fc,)) )) ’ 

where K, II and E are the complete elliptic integrals 123 § 19.2(h)], while E is the unperturbed highest-energy state 
energy satisfying the quantization condition, 

^ - tt(1 -z + )-t{E< A/2) = (l - x) ^ ( K(fc) - n(a2 ’ k) ) “ ~~ +E(fc) ’ (10) 


with 1 (•) denoting the indicator function. 

These complicated expressions constitute a solution to 
the problem of semiclassical quantization but offer little 
insight into the dimer’s behavior. Nonetheless, some of 
the problem’s structure has become apparent: 

1. The splitting depends on E and A only through the 
turning points z± and the combination (1 — 2E/A). 
The sign of this last quantity distinguishes between 
the two types of motion depicted in Figure [6j 1 — 
2E/A > 0 for rotations (orbits surrounding one 
of the poles at 2 = ±1) and 1 — 2E/K < 0 for 
librations. 

2. The only nonelementary functions in the expres¬ 
sions above are the complete elliptic integrals K, 
E, and II. When they do appear they all take the 
same argument (modulus), either k or k ', which is 
a measure of the size of the classical orbit. 


This structure can be exploited to find much simpler ex¬ 
pressions for the splitting, valid in the limit of A > 1. 

Let us first rescale the energy through a linear trans¬ 
formation: 



(A-l ) 2 

2A 


( 11 ) 


The rescaled energy e lies in [0, 1) for any orbit in the self¬ 
trapping region. The highest-energy state orbit has an 
area h/ 2, while the total semiclassical action of a dimer 
with N particles is h(N + 1). As N increases, both the 
highest-energy state energy e and the dimensionless mea¬ 
sure of orbit size k [Eq. ([8])] become small. If the highest- 
energy state orbit is a libration (e < (A — l) -2 ), expand¬ 
ing Eq. (101 to lowest order in k and e and solving for e 
gives an estimate of the highest-energy state energy, 


2AVA 2 - 1 
(A — 1) 2 (N + 1)' 


(12) 
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FIG. 7. Comparison of semiclassical estimates of the splitting 
with exact diagonalization. The analytical approximation of 
Eq. (131 (red line) agrees closely with the results of exact 
diagonalization (blue dashed line). In contrast, the approxi¬ 
mation of [20] (green dot-dashed line) performs poorly in this 
low-A regime, especially for larger N. The black vertical line 
marks the A value below which the semiclassical approxima¬ 
tion breaks down because the area of phase space associated 
with the self-trapped region is less than h/2. 


This estimate is very good: the relative error in approxi¬ 
mating the numerical semiclassical result is less than 1% 
for N = 20 and A = 1.25, and decreases with both N 
and A. Analogous expansions for the classical orbital 
frequency and the tunneling phase lead to the following 
expression for the ground state splitting: 


A E 




LO 


(JV+lXt—e) 


(13) 


where zq 


1 - 


l 

A 2 


is the position of the self-trapping 


fixed point and co = VA 2 — 1 is the frequency of motion 
about it. The tunneling frequency A E/h decreases ex¬ 
ponentially with the “barrier width” ss zq, the “barrier 
height” ss (1—e) and the number of atoms N. The details 
of the calculation are described in Appendix [C| 

Figure [7] compares the semiclassical splitting estimates 
with the results of exact diagonalization of the Bose- 
Hubbard model. The results of solving the quantiza¬ 
tion problem numerically are not shown: except for A 


so small that not even one semiclassical orbit fits within 
the self-trapping region, they agree very closely with the 
exact Bose-Hubbard splitting. The analytic approxima¬ 
tion discussed in this section is generally within a factor 
of 2 of the exact result, and improves with N. Since the 
splitting changes by as many as 15 orders of magnitude 
over the investigated range of A, this agreement amounts 
to remarkably robust performance. 

A different closed-form semiclassical approximation to 
A E was obtained in m and refined in £20]. This last ap¬ 
proximation attains an excellent accuracy, on the order 
of a few percent, but only for U ~ J. In the context of 
cold atomic experiments, in which the atom number is on 
the order of hundreds, this corresponds to astronomically 
small tunneling frequencies (well below 10“ 100 Hz). For 
U <C J, where the tunneling frequency becomes large, the 
approximation of [201 is many orders of magnitude from 
the true value (see Figure [7]). Therefore, the approxi¬ 
mation we provide in Eq. (13) is the first closed-form 
expression valid in the experimentally relevant regime. 


IV. DISCUSSION 

In this section, we consider the implications of the 
analysis presented above for three problems: determining 
the time scale for macroscopic entanglement, producing 
quantum speedup of dissipation, and obtaining experi¬ 
mental confirmation. 


A. Time scale for macroscopic entanglement 

Tunneling in the self-trapping regime leads to the 
generation of entangled superpositions of many-particle 
states, or macroscopic entanglement m • The entangle¬ 
ment between the two modes is maximized at times T /4 
and 3T/4, where T is the tunneling period. Therefore, 
our semiclassical estimate of the tunneling frequency im¬ 
mediately yields an estimate of the time required for en¬ 
tanglement generation. It is notable that the dynamics 
of entanglement, a profoundly unclassical phenomenon, 
is captured by the first quantum correction to the (clas¬ 
sical) BJJ model. 


B. Quantum speedup of dissipation 

So far we have considered only an isolated Bose- 
Hubbard dimer. We will now discuss the central role 
tunneling in the self-trapped regime plays in a quantita¬ 
tively accurate model of a dissipative dimer that leaks 
atoms to the environment. 

Consider a coherent state of N bosons centered at 
one of the self-trapping fixed points, say the left well. 
We will attempt to model its dynamics within a two- 
dimensional subspace of the full system’s Hilbert space, 
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the subspace spanned by the symmetric and antisym¬ 
metric energy eigenstates, \E$) and | Ea)- In the basis of 
states localized in the two wells, | 1 ) = (\Es) + \Ea))/V 2 
and | 2 ) = (| E$) — \Ea))/V2, the Hamiltonian is repre¬ 
sented by the matrix, 

f E A E\ 

\AE E ) 


where E = (E$ + Ea )/2 and A E = (Es — Ea)/2. 
These parameters can be calculated semiclassically with 
high accuracy as we have shown in the preceding section 
[Eq. (12 1 and Eq. (13)], though we use exact values in 
the simulation discussed below. The initial condition is 
the localized state |1). Now, assume there is decay from 
the right well at a rate 7 . In the two-level model this is 
described by the effective decay rates, 


Ti = —7 (11 at 


A I A 


leading to the effective Hamiltonian, 

(2) _(E- iTi/2 A E \ 

n oS - y A E E- iT 2 / 2 ) ' 


(14) 


This simple model can be used to estimate how the prob¬ 
ability of all N atoms remaining in the system diminishes 
over time. To evaluate the results, we compare them to 
those obtained using the complete coherent state and the 
full master equation |38l - t40] . 

P = ~t[H, p] - ^ (a\aip + pa\ai - 2a 1 pa[ S j . (15) 

The probabilities of remaining in the N atom subspace 
predicted using the two Hamiltonians are shown in Fig- 
ure[ 8 j If many-body tunneling between the fixed points is 
neglected (A E = 0 ), the rate of atom loss is significantly 
underestimated. But when the correct value of A E is 
used, the effective two-state model produces results al¬ 
most indistinguishable from the full Bose-Hubbard. Re¬ 
markably, we can thus reproduce the decay dynamics of a 
correlated many-body system using only two parameters, 
E and A E, which can be calculated semiclassically. 


C. Prospects of experimental observation 



FIG. 8. Correctly estimating the rate of tunneling between 
the self-trapping fixed points is critical to predicting the atom 
loss rate from a leaky dimer. The probability of finding all 
N atoms in the system over time is plotted for three different 
mod els. The dashed green line is the simple Hamiltonian of 
Eq. (141, based only 011 two parameters E and AE which can 
be calculated semiclassically. It overlaps with the numerically 
exact results obtained by integrating the many-body master 
equation of Eq. ( |15| ) (solid blue line). The simple model with 
A E set to zero differs significantly (dotted red line). (J = 
1 Hz, [7 = 4/5 Hz, N = 6) 



FIG. 9. Frequency of tunneling between the fixed points ver¬ 
sus A for U = 2n x 0.063 Hz and N = 500, the experimental 
parameters of m The mean-field prediction is also shown 
for reference. 


The BJJ dynamics of the BEC dimer was experimen¬ 
tally mapped out in great detail a few years ago HU. 
Could a similar experiment observe tunneling between 
the fixed points for A > 1 ? 

Experimental realizations of the dimer fall into two 
categories: “external” and “internal” m, or those uti¬ 
lizing two spatially separated wells and those using two 
internal states of atoms. Correctly describing the dynam¬ 
ics of the spatially separated wells requires going beyond 
the Bose-Hubbard model that was our starting point in 
this work, as the localized orbitals associated with the 
operators a t . a\ are time-dependent |42| . Fortunately, 


this complication does not arise in the case of internal 
states [43]. Therefore, the tunneling and dissipation en¬ 
hancement effects we have described are most likely to 
be observed in experiments relying on internal states. 

The expected tunneling frequency given the experi¬ 
mental parameters of m is shown in Figure [9] The 
frequency is on the order of a few Hertz. Since the atom 
decay times reported in this experiment are ~ 100 ms, 
the tunneling frequency is too small to be observed at 
present. However, an order of magnitude improvement 
in atom retention times would render experimental ob¬ 
servation feasible. 
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At first glance, it may seem that the retention time 
limitation could be sidestepped by lowering both N and 
J by the same factor. Since the quantum tunneling 
time depends on N exponentially, but on J only linearly 
[Eq. @] , this could speed up the semiclassical dynamics 
while keeping A constant. Unfortunately, the experiment 
of El was already carried out at the lowest J currently 
accessible: lowering it even more introduces unacceptable 
noise due to EM fluctuations [33]. 
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FIG. 10. Graphical representation of the roots of Eq. (A11. 


D. Beyond the dimer: semiclassical quantization 
for lattices 

Although our analysis was limited to the dimer, anal¬ 
ogous processes should occur in a system with multiple 
states, only one of which has an appreciable population. 
The Bose-Hubbard Hamiltonian can be straightforwardly 
extended to such systems; in the case of the trimer, self¬ 
trapping has been demonstrated in both the quantum 
model and its classical limit [45], 46]. However, carry¬ 
ing out semiclassical quantization is difficult because the 
classical model is now chaotic. So far, progress has only 
been made for the case of very small and very large 
J/U [47] , i.e. precisely the region of parameter space 
where tunneling between the self-trapping points does 
not take place. Therefore, the extension of our results 
beyond the dimer is likely to prove challenging. 


V. SUMMARY & OUTLOOK 

We have studied the tunneling between the self- 
trapped fixed points of the BEC dimer using a semiclas¬ 
sical approach. We derived an exact solution to the prob¬ 
lem in terms of elliptic integrals giving the phase space ar¬ 
eas of semiclassical orbits. For particle numbers N 1, 
the semiclassical ground state orbit and (appropriately 
transformed) energy become small; in this limit we found 
an approximate closed-form expression for the tunneling 
frequency that is accurate in the experimentally relevant 
parameter range. The tunneling frequency decreases ex¬ 
ponentially with the effective width and height of bar¬ 
riers in phase space, as well as the number of particles. 
Nonetheless, accounting for the tunneling is crucial to 
obtaining quantitatively accurate estimates of atom loss 
rates in a leaky dimer. 
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Appendix A: Derivation of Equation (|7j) 


In this appendix, we use Eq. ([5]) , the quantization con¬ 
dition of Graefe and Korsch [l7l . to derive an approx¬ 
imate expression for the energy splitting of the nearly- 
degenerate self-trapped eigenstates. This expression and 
its derivation have been known to scholars of the WKB 
approximation (see |1%] . p. 49, or [3B], p. 52), but the 
discussion we give here is more complete than that found 
in other sources. 

Eq. <© can be rewritten as, 


cos ( 2S W S ^) 


1 

a/ 1 + exp(27rS' £ ) 


(Al) 


Considered as a function of x = 2 S w — S^, this equa¬ 
tion has pairs of solutions symmetrically spaced about 
(2 n+ l)7r (see Figure [To|. The pairs of roots coalesce as 
S e —>■ —00: in the absence of tunneling, states come in 
degenerate pairs, one localized in each well. Let the two 
solutions near x = 7r be x±, with x+ > n and X- < 7 r. 
We have, 


t = = 

COS x± 

where the sign difference on the right-hand-side arises 
because sin(x) changes sign at x = 7r, between x __ and 
x + . 

Recall that x = 2 S w — is a function of energy. As¬ 
sume the ground state energy splitting A E is sufficiently 
small that x{E) is approximately linear in an interval of 
width A E about the ground state energy, E 0 . Then, 


= x{Eq ± AE/2), 


and Eq. (A2) gives, 


tan(2 S W (E 0 ± AE/2) - S 0 (S e (£b ± AE/2))) = 

=F exp(7rS' £ (E 0 ± AE/2)), 


or, 

2S W (E 0 ± AE/2) - S^(S £ (Eq ± AE/2)) = 

=F arctan (exp nS e (Eo ± AE/2)). 













Expanding to first order about Eq, 


Appendix B: Elliptic integral expressions for T, S w 
and S f 


2 S,„ - ± 2 


.8S„ 


dS 4 , 
dS e 


dS e 

~dE 


dE 

+ arctan(exp(7rS' e )) — 


A E _ 
2 “ 
27T 


dS e A E 


cosh(7rS' e ) dE 2 

Subtracting the lower signs from the upper signs and re¬ 
arranging yields, 

A E arctan exp (7^) 

~T~ “ ~ 


J dE 


ds^ es„ 

dS e dE 


(A3) 


Consider the second term in the denominator. Letting 
£ = S e and using the definition of [Eq. ([6])], the unit¬ 
less derivative can be written as. 


dS 6 


1 


1 


1 


where ip is the digamma function, defined as 

r'C t) 


(A4) 


xp{t) = 


m 


For \t\ > 3, excellent approximation (good to 0.03%) to 
this function is provided by the asymptotic expansion 
5.11.2], 


ip(t) ? 

Using this expansion, 


In t — - - 

2 t 


12 1 2 ' 


dS ± 
dS. 


1 


In 1 + 


4^ 

3 — 8£ 2 (1 - 2£ 2 ) 
24£ 2 (1 + 2£ 2 ) 2 


4 1 + 2£ 2 
3(1+ 4£ 2 ) 2 


This expression is already smaller than 0.01 at £ = 2, and 
decreases with £ as l/£ 2 . Since the phase space deriva¬ 
tives dS w jdE and dS e /dE are of the same order, and 
£ = S e is of order N, the second term in the denomina¬ 
tor of Eq. (A3) can be neglected: 

arctan exp irS e 


A E = -- 


as+ 

dE 


Since the splitting is small, exp(7r5' e ) <C 1 and so 
arctan exp (jrS e ) ~ exp(7rS' e ). If we let T = 2tt/uj be 
the period of the orbit corresponding to the action 2 S w , 


0 dS w 

dE 

Neglecting the second 
Eq. (A3), we get, 


_ 1 _ 2ti- 

h huj 

term in the denominator of 


hco 


A E = -exp(7r5' e ). 


(A5) 


The negative sign of A E indicates that x + is actually 
lower in energy than X-. 

As a special case, this result applies to a single particle 
in a double-well potential described by the Schodinger 
equation. For that special case there exist a simpler 
derivation of Eq. (A5): see [JSj, §50. 


To perform actual calculations using the formula, 

A E = — exp(7r S e ), 

TV 

we need to find explicit expressions for to (or the corre¬ 
sponding period T) and S € in terms of E and A. It will 
also prove useful to find an expression for 2 S w , the ac¬ 
tion associated with the self-trapped orbit, which deter¬ 
mines the energy about which the splitting takes place. 
All of these quantities depend on the shape of the classi¬ 
cal orbits of the mean-field Hamiltonian of Eq. ([3]). The 
equation of the orbit is, 

Az 2 — 2 E 

(p{z, E, A) = arccos ^ ^ , (Bl) 


and the classical turning points of the orbits (see Fig¬ 
ure [6]) are, 


z±(E,A) = 


I ±Vl~2EA + A 2 + AE-1 

AV2 ' 


(B2) 


In what follows, we will generally suppress the explicit 
dependence of <j> and z± on E and A to obtain clearer 
expressions. Recall that we defined the dimensionless 
measure of orbit size as, 



We begin with the simplest problem, that of deriving 
an expression for the orbit period T. The approach to 
computing the action integrals S e and S w is the same, 
but the technical details are more involved. 

See [50] and the references therein for a deeper look at 
the geometry of the classical model and its relationship 
to Bose-Hubbard dynamics. 


1. Period T of the classical orbit 


The equation of motion for z is, 


dH n— 


and so the period is, 


T = 2 


r z + dt , 

— dz 
, dz 


r z + 


= 2 


dz 


z_ \J\ — z 2 sin (p(z) 


(B3) 


(B4) 


Since sin(arccosa;) = y/l — x 2 , we can use Eq. (Bl) to 
eliminate the trigonometric functions: 


T = 4 


r 

J Z- 


dz 


\J 4(1 - z 2 ) - (Az 2 - 2 E ) 2 ' 


(B5) 
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Although at first glance this expression has a very com¬ 
plicated structure, the polynomial in the denominator 
(which is also encountered in the S w and S e integrals) 
can be rewritten in the more suggestive form, 


4(1 - z 2 ) - {Az 2 - 2 E) 2 = -A 2 {z 2 - zl)(z 2 - z 2 ). 
The period is therefore, 


4 r z + dz 

^' z - (z* - z 2 ){-zl + z*) 



4 

A z+ 


K (k), 


(B6) 


(B7) 


where K is the complete elliptic integral of the first kind. 
Note that in this expression, time is measured in the 
dimensionless units introduced with the Hamiltonian of 
Eq. (§. Converting the units to seconds, 


T = - 


J Az. 


-K(fc), 


(B8) 


where J is measured in Hertz. 


2. Action of the classical orbit 

The phase space areas (and so actions) associated with 
the classical orbits can be found by integrating </>(z). For 
an orbit in the self-trapping region, the action is 


S(E, A) = h 


N + l 

47r 


^2 J 7r — cf>(z) dz 


(B9) 


+ 2tt(1 - z+) 1 (E < A/2) 


The prefactor h^AA normalizes the total area of phase 
space to be h(N + 1), with N the number of particles. If 
E < A/2, the orbit is a rotation orbit (see Figure [6]) and 
the area of the “cap” at |z| > z + is added to the integral 
of <j)(z). 


The integral in Eq. (B91 can be simplified through an 
integration by parts: 

z 2 (z 2 + 2^A) 


r z + r z + 

/ 7r — (f>{z ) dz= — 

Jz - J z - (l-z^^(z^z 2 )(-z^ + z 2 + ) 

where the boundary term is zero since < f>(z±) = ir/2. This 
is an elliptic integral E3 § 19.2(i)] and can be reduced to 
the canonical elliptic integrals using a partial fraction 
decomposition. Let, 


P = ~(* 2 - z\){z 2 - z 2 _). 


Then, 


/ ■ z + / 2>E\ 1 

7T — (f>(z) dz = — 2L|-E(fc) + f 1 ——J —X 


K(k)- 


1-4 


n (a 2 , fc) 


(BIO) 


where K (k), E(fc), and n(a 2 , k) are complete elliptic inte¬ 
grals of the first, second and third kinds, k is the measure 
of orbit size defined in Eq. ([8]), and 


a 


2 



3. Tunneling action S e 


The “tunneling action” is defined analogously to the 
orbit action, 

N+l f z -( E A) 

Se(E, A) =- —— — • 2 / |tt — <j>(z, E, A)| dz, 

4?r A) 

with the absolute value necessary because <f>(z , E, A) may 
be complex within the region of integration. In fact, in 
the self-trapping region (E > 1, A > 1) the argument 
of the arccosine in <f>(z,E, A) is smaller than —1 for all 
2 € [-Z-, z_]. Consequently, taking advantage of the 
identity, 

arccos(—1 — x) = n — * arccosh (1 + a;), 


one may rewrite S e as, 


= - 


N+l 


arccosh 


2 E - Az 2 \ 

2 VT 77 ^ 2 J 


dz. 


As in the case of the orbit action, S e can be recast as an 
elliptic integral through integration by parts, and then 
reduced to a sum of canonical elliptic integrals using a 
partial fractions expansion. The result is, 


TrS e 
N + 1 


= -( 1 -x)£ n( *A*' ) 

+ ^ + (K(k')-E(k')), 


(BH) 


where k 1 = \J 1 — k 2 and we have used identity 19.6.5 
in |3Z| . 


Appendix C: Approximate solution to the 
quantization problem for large N 

In this section, we derive an approximate semiclassical 
expression for the splitting by expanding the integrals of 
the previous section in small orbit sizes, k, and energies, 

e. 


1. Approximate orbit frequency 

To lowest order, 

< ci > 

A higher-order expansion is unnecessary because A E de¬ 
pends on e primarily through the tunneling phase in the 
exponent. 























10 


2. Energy of the highest-energy state 


error 

0.008 r 


Many of the quantities encountered in our discussion 
so far can be expressed more simply in terms of e [the 
normalized energy relative to the maximum of E —see 
Eq. (Ill] than E. For instance, the classical turning 
points are, 

z± = l-^(lT(A-l)v^) 2 

and the dimensionless measure of orbit size is, 

- A 4^ 


r = 


(y/e + l) 2 + A(1 - e)' 


The quantization condition of Eq. (101 reads 


— 7r(l — z+) • 1 (e > (A — 1) 2 ) = — z+E(k) 

1 


N + l 

1 - (A - l) 2 e 1 


A 2 


z+ 


K(jfe) - 


l- z\ 


n(a , k) I , (C2) 


with k and z± given by the expressions in the previ¬ 
ous section. Consider the case e < (A — l) -2 , when the 
highest-energy state orbit is a libration. Expanding the 
elliptic integrals to lowest order in k and then to lowest 
order in e m and then solving for e gives a first-order 
estimate of the highest-energy state energy, 


e = 


2AVA 2 - 1 
(A - l) 2 (N + 1)' 


(C3) 


As was already remarked in the main text, this estimate 
is very good. See also Figure [TT| 

What happens if the nonlinearity is sufficiently high 
that the highest-energy state orbit is a rotation (i.e., 
(A — 1)~ 2 < e <C 1)? It turns out that this case cannot 
be successfully treated using the same approach. The 
term t _ z2 n(a 2 , k) becomes ill-behaved, with both the 

prefactor and a 2 very large. The terms of the small-fc 2 
expansion of II(a 2 ,A:) are proportional to powers of a 2 
[cf. Eq. 19.5.4 in [37]], so keeping only the lowest-order 
terms in k 2 is no longer legitimate. But the difficulty of 
extending our semiclassical method to this part of the 
parameter space is not a major concern, for two reasons: 

1. The nonlinearity required for the ground state orbit 
to enclose the point z = 1 is large indeed, especially 
for larger atom numbers. From Eq. (12), the con¬ 


dition e > (A — l) 2 can be estimated to imply, 


2A\/A 2 - 1 « 2A 2 > N + 1. 


(C4) 


2. The limit of very strong nonlinearity is particu¬ 
larly easy to treat using quantum perturbation the¬ 
ory [HI [301 [331 [34] . 


0.006 

0.004 

0.002 


A 


0 20 40 60 80 100 

error 
0.0501 • 

0.010 * 

0.005 


0.001 

5x 10“ 4 


N 


20 40 60 80 100 

FIG. 11. Relative error in approximating the (numerically ex¬ 
act) solution of Eq. ( |C2[ ) with the lowest-order approximation 
of Eq. (C31. The figure on the left shows the dependence on 


A (for N = 20) and that on the right—the dependence on N 
(for A = 2). 


3. Approximate tunneling action 

Finding a good large-A approximation for the tun¬ 
neling action [Eq. (Bill] is more difficult because both 


n(z + 2 ,fc') and K (k') — E (k r ) diverge in the limit k' = 
y/l — k 2 -A 1_. The lowest order asymptotic approxima¬ 
tion is of O(e 0 ): 


7 tS £ 

' N + 1 


y\T~ 


+ In (A + VA 2 - l) ■ 


It is possible to derive higher-order approximations by 
combining the known asymptotic expansions of the com¬ 
plete elliptic integrals, but they are complex and disap¬ 
pointingly inaccurate, except for large N and either very 
large or very small A. 

Instead of pursuing a formal expansion, let’s attempt 
an ad hoc improvement of the zeroth-order expression. 
S € is a measure of the barrier to tunneling; as the ground 
state approaches the separatrix (e -A 1), the barrier 
should disappear. The simplest way to enforce this be¬ 
havior is to multiply the O(e 0 ) expression by (1 — e): 


7 tS € 
N+ 1 


Va^T 

A 


+ In (A+ VA 2 -l)^J (1 


-e). 

(C5) 






























bound on the magnitude of S e for all A. 
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Log 10 AE 



Log 10 AE 



FIG. 12. Approximations to the semiclassical highest-energy 
state splitting. Numerical solutions to Eq. © are shows as 
blue dots; Eq. (C 61 is plotted as the solid red line, while 
Eq. ‘ 


( |C6| ) with e = 0 is shown in dashed blue. The black 
vertical line marks the point where the semiclassical approx¬ 
imation must break down because the area of phase space 
associated with the self-trapped region is less than h/2. 


This ansatz works remarkably well; furthermore, unlike 
the asymptotic expansions which may be either smaller 
or larger than the true value, Eq. (C5) gives an upper 


4. Approximate splitting formula 

By combining the approximate expressions for the clas¬ 
sical orbital frequency and the tunneling phase, we arrive 
at the following expression for the highest-energy state 
splitting: 


AE 


Hco 

7T 



UJ 


(JV+l)(l-e) 


where Zq = ^Jl — is the position of the classical po¬ 
tential maximum and u> = aJ 1 — is the frequency of 
motion about it [cf. Eq. (Cl)]. In this expression, the 
frequency is measured in the dimensionless units intro¬ 
duced with the Hamiltonian of Eq. ^. In the units of J 
and U (Hertz), 


T w / 1 

7T \ U! 


AE « 2 J- ( —e~ z ° ) 


yw+l)(l—e) 


(C6) 


Figure [12] shows a comparison of this approximation 
with the numerical solution of the semiclassical quanti¬ 
zation condition [Eq. ([5|]. Since our approximation to S e 
overestimates the barrier to tunneling, the tunneling fre¬ 
quency is generally underestimated, except close to the 
bifurcation where the dependence of to on e (which we 
neglect) becomes important. Some qualitative features 
of the dependence of AE on A can be reproduced even 
without the factor of (1 — e) in the exponent, and the 
agreement with the numerical solution improves as N in¬ 
creases. However, this e = 0 approximation to AE is 
generally not within an order of magnitude of the nu¬ 
merically computed value. 
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